#include "FourierTron.h"
#include <iostream>



using namespace std;
using namespace SeamStress;


typedef unsigned long int ulong;


FourierTron::FourierTron(vector<Seamstress*> *ss, unsigned long int nt, mpfr_prec_t p, unsigned long int cache) : pins(this, ss), cache_size(cache)
{
  if(nt==0){nthreads=1;}
  else{nthreads=nt;}
  if(ss==NULL){nthreads=1;}
  else{if(nthreads > ss->size()){nthreads = ss->size();}}
  
  prec = p;
  
  
  pi.set(0,0,prec);
  mpfr_const_pi(pi.real_part, GMP_RNDN);
  tempval1.set(0,0,prec);
  tempval2.set(0,0,prec);
  
  
  threadvec1.assign(nthreads, tempval1);
  threadvec2.assign(nthreads, tempval1);
  
  localvec.assign(nthreads, tempvec1);
}


FourierTron::~FourierTron()
{
  
}



void FourierTron::inverseCT_FFTThread1(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  for(unsigned long int i=w;i<thread_var1;i+=thread_tot)
  {
    (*tempvecpointer1)[i].mul(tempval1);
  }
}


void FourierTron::inverseCT_FFT(vector<cxfloat> &input, unsigned long int N)
{
  for(unsigned long int i=0;i<N;i++){input[i].conj();}
  CT_FFT(input, N);
  for(unsigned long int i=0;i<N;i++){input[i].conj();}
  tempval1.set(1,0);
  tempval2.set(N,0);
  tempval1.div(tempval2);
  thread_var1=N;
  if(N<nthreads){thread_tot=N;}
  else{thread_tot=nthreads;}
  tempvecpointer1 = &input;
  
  pins.sewStraight(&FourierTron::inverseCT_FFTThread1, thread_tot);
}


void FourierTron::inverseCT_FFT(const vector<cxfloat> &input, vector<cxfloat> &output, unsigned long int N)
{
  for(unsigned long int i=0;i<N;i++){output[i].set(input[i]);}
  inverseCT_FFT(output, N);
}


void FourierTron::CT_FFT_mul3threads(cxfloat_vec& input, unsigned long int N)
{
  CT_FFT_mul3threadsSplit_DATA.input = &input;
  CT_FFT_mul3threadsSplit_DATA.N = N;
  cxfloat_vec vec0;
  cxfloat_vec vec1;
  cxfloat_vec vec3;
  cxfloat_vec vec4;
  
  CT_FFT_mul3threadsSplit_DATA.vec0 = &vec0;
  CT_FFT_mul3threadsSplit_DATA.vec1 = &vec1;
  CT_FFT_mul3threadsSplit_DATA.vec3 = &vec3;
  CT_FFT_mul3threadsSplit_DATA.vec4 = &vec4;
  unsigned int tthreads = 4;if(nthreads < 4){tthreads = nthreads;}
  CT_FFT_mul3threadsSplit_DATA.tthreads = tthreads;
  pins.sewStraight( &FourierTron::CT_FFT_mul3threadsSplit, tthreads );
  {
    vector<SeamStress::Seamstress*>* ss = pins.seamstresses;
    vector<SeamStress::Seamstress*> one_ss_1;
    one_ss_1.push_back( (*ss)[0] );
    SeamStress::Pincushion<FourierTron> one_pin_1(this, &one_ss_1);
    vector<SeamStress::Seamstress*> third_1_ss;
    for(ulong i=1, loopsizei=(2*(ss->size())/3);i<loopsizei;i+=1){third_1_ss.push_back( (*ss)[i] );}
    third_1_ss.push_back(nullptr);
    CT_FFT_subFFT_TYPE CT_FFT_subFFT_DATA_1;
    CT_FFT_subFFT_DATA_1.ss = &third_1_ss;
    CT_FFT_subFFT_DATA_1.input = &vec0;
    vector<void*> threadvec_1;threadvec_1.push_back( (void*)(&CT_FFT_subFFT_DATA_1) );
    one_pin_1.sewOpenly( &FourierTron::CT_FFT_subFFT, threadvec_1 );
    
    vector<SeamStress::Seamstress*> third_3_ss;
    for(ulong i=(2*(ss->size())/3), loopsizei=(ss->size());i<loopsizei;i+=1){third_3_ss.push_back( (*ss)[i] );}
    CT_FFT_subFFT_TYPE CT_FFT_subFFT_DATA_3;
    CT_FFT_subFFT_DATA_3.ss = &third_3_ss;
    CT_FFT_subFFT_DATA_3.input = &vec1;
    CT_FFT_subFFT((void*)(&CT_FFT_subFFT_DATA_3));
    
    one_pin_1.tieOff();
  }
  CT_FFT(vec3, vec3.size());
  
  CT_FFT_subTwiddle_DATA.even = &vec0;
  CT_FFT_subTwiddle_DATA.even1 = &vec1;
  CT_FFT_subTwiddle_DATA.odd1 = &vec3;
  CT_FFT_subTwiddle_DATA.output = &vec4;
  pins.sewStraight(&FourierTron::CT_FFT_subArrange, nthreads);
  CT_FFT_subArrange_DATA.input = &vec4;
  CT_FFT_subArrange_DATA.output = &input;
  CT_FFT_subArrange_DATA.N = N;
  pins.sewStraight(&FourierTron::CT_FFT_copyArrange, nthreads);
  
  CT_FFT_subTwiddle_DATA.output = &input;
  CT_FFT_subTwiddle_DATA.N = N;
  pins.sewStraight(&FourierTron::CT_FFT_subTwiddle, nthreads);
}


void FourierTron::CT_FFT_copyArrange(void* arg)
{
  ulong w = *((ulong*)arg);
  
  cxfloat_vec& input = *(CT_FFT_subArrange_DATA.input);
  cxfloat_vec& output = *(CT_FFT_subArrange_DATA.output);
  unsigned long int N = CT_FFT_subArrange_DATA.N;
  
  for(ulong i=w,loopsizei=N;i<loopsizei;i+=nthreads)
  {
    output[i].point(input[i]);
    input[i].set(nullptr);
  }
}


void FourierTron::CT_FFT_subArrange(void* arg)
{
  ulong w = *((ulong*)arg);
  
  cxfloat_vec& even = *(CT_FFT_subTwiddle_DATA.even);
  cxfloat_vec& odd1 = *(CT_FFT_subTwiddle_DATA.odd1);
  cxfloat_vec& even1 = *(CT_FFT_subTwiddle_DATA.even1);
  cxfloat_vec& output = *(CT_FFT_subTwiddle_DATA.output);
  
  for(ulong i=w,loopsizei=(output.size()/4);i<loopsizei;i+=nthreads)
  {
    output[i].point( even[i] );
    output[i + loopsizei ].point( even[i + loopsizei] );
    output[i + 2*loopsizei].point( even1[i] );
    output[i + 3*loopsizei].point( odd1[i] );
    
    even[i].set(nullptr);
    even[i + loopsizei].set(nullptr);
    even1[i].set(nullptr);
    odd1[i].set(nullptr);
  }
}


void FourierTron::CT_FFT_subTwiddle(void* arg)
{
  ulong w = *((ulong*)arg);
  
  cxfloat_vec& output = *(CT_FFT_subTwiddle_DATA.output);
  unsigned long int N = CT_FFT_subTwiddle_DATA.N;
  
  ulong log2N = (ulong)(__builtin_ctzl(N)) - 1;
  
  cxfloat twidval(0,0,prec);
  cxfloat newodd1(0,0,prec);
  cxfloat newodd2(0,0,prec);
  cxfloat E(0,0,prec);
  cxfloat EE(0,0,prec);
  cxfloat E1(0,0,prec);
  cxfloat O1(0,0,prec);
  for(ulong i=w, loopsizei=(N/4);i<loopsizei;i+=nthreads)
  {
    cxfloat& X1 = output[i];
    cxfloat& X2 = output[i + loopsizei  ];
    cxfloat& X3 = output[i + 2*loopsizei];
    cxfloat& X4 = output[i + 3*loopsizei];
    E.set(X1);
    EE.set(X2);
    E1.set(X3);
    O1.set(X4);
    
    twiddleVal(log2N-1, i, twidval);
    newodd1.set(O1);newodd1.mul(twidval);
    newodd2.set(newodd1);newodd2.neg();
    newodd1.add(E1);newodd2.add(E1);
    twiddleVal(log2N, i, twidval);
    X1.set(newodd1);X1.mul(twidval);
    twiddleVal(log2N, i+loopsizei, twidval);
    X2.set(newodd2);X2.mul(twidval);
    X3.set(X1);X3.neg();
    X4.set(X2);X4.neg();
    X1.add(E);X2.add(EE);X3.add(E);X4.add(EE);
  }
}


void FourierTron::CT_FFT_subFFT(void* arg)
{
  CT_FFT_subFFT_TYPE& data = *((CT_FFT_subFFT_TYPE*)arg);
  
  FourierTron ftron(data.ss, data.ss->size(), prec, cache_size);
  ftron.CT_FFT(*(data.input), data.input->size());
}


void FourierTron::CT_FFT_mul3threadsSplit(void* arg)
{
  ulong w = *((ulong*)arg);
  
  const cxfloat_vec& input = *(CT_FFT_mul3threadsSplit_DATA.input);
  unsigned long int N = CT_FFT_mul3threadsSplit_DATA.N;
  cxfloat_vec& vec0 = *(CT_FFT_mul3threadsSplit_DATA.vec0);
  cxfloat_vec& vec1 = *(CT_FFT_mul3threadsSplit_DATA.vec1);
  cxfloat_vec& vec3 = *(CT_FFT_mul3threadsSplit_DATA.vec3);
  cxfloat_vec& vec4 = *(CT_FFT_mul3threadsSplit_DATA.vec4);
  unsigned int tthreads = CT_FFT_mul3threadsSplit_DATA.tthreads;
  
  for(ulong i=w, loopsizei=N;i<loopsizei;i+=tthreads)
  {
    if(i == 0)
    {
      vec0.reserve( N/2 );
      for(ulong j=0, loopsizej=N/2;j<loopsizej;j+=1)
      {
        vec0.emplace_back(nullptr);
        vec0.back().point( input[2*j] );
      }
    }
    if(i == 1)
    {
      vec1.reserve( N/4 );
      for(ulong j=0, loopsizej=N/4;j<loopsizej;j+=1)
      {
        vec1.emplace_back(nullptr);
        vec1.back().point( input[4*j+1] );
      }
    }
    if(i == 2)
    {
      vec3.reserve( N/4 );
      for(ulong j=0, loopsizej=N/4;j<loopsizej;j+=1)
      {
        vec3.emplace_back(nullptr);
        vec3.back().point( input[4*j+3] );
      }
    }
    if(i == 3)
    {
      vec4.reserve( N );
      for(ulong j=0, loopsizej=N;j<loopsizej;j+=1)
      {
        vec4.emplace_back(nullptr);
      }
    }
  }
}


void FourierTron::CT_FFT(vector<cxfloat> &input, unsigned long int N)
{
  unsigned long int chunk_size = cache_size/(prec>>1);
  step_depth = 0;
  while(chunk_size != 0){chunk_size>>=1;step_depth++;}
  //check that N is a power of two
  try{if((N & (N-1)) != 0){throw (char*)("second argument must be a power of two");}}
  catch(char* str){cout<<"Exception from FourierTron::CT_FFT : "<<str<<endl;throw;return;}
  //check that input.size() >= N
  try{if(input.size() < N){throw (char*)("input vector size smaller than second argument");}}
  catch(char* str){cout<<"Exception from FourierTron::CT_FFT : "<<str<<endl;throw;return;}
  
  if( ((nthreads & (nthreads-1)) != 0) )
  {
    if( (nthreads%3 == 0) && (N >= 32) && ( N > 16*nthreads ) )
    {
      CT_FFT_mul3threads(input, N);
      return;
    }
  }
  
  fft_size=N;
  calculateTwiddle(N);
  unsigned long int logN=0;
  unsigned long int NN = (N>>1);
  while(NN!=0){logN++;NN = (NN>>1);}
  bitReverse(logN);
  tempvecpointer1 = &input;
  thread_tot=1;
  NN=0;
  while(thread_tot<nthreads){thread_tot<<=1;NN++;}
  if(thread_tot>nthreads){thread_tot>>=1;NN--;}
  N>>=1;
  while(thread_tot>N){thread_tot>>=1;NN--;}
  thread_var1 = logN;
  pins.sewStraight(&FourierTron::CT_FFT_swap, thread_tot);
  thread_var2 = thread_tot;
  NN = fft_size/thread_tot;
  while(NN<thread_tot)
  {
    NN<<=1;
    thread_var2>>=1;
  }
  NN = thread_var2>>1;
  thread_var3=0;
  while(NN!=0){NN>>=1;thread_var3++;}
  if(step_depth==0){pins.sewStraight(&FourierTron::CT_FFT_step, thread_tot);}
  else{pins.sewStraight(&FourierTron::CT_FFT_step_cache, thread_tot);}
  while(thread_var2<thread_tot)
  {
    thread_var4 = thread_tot;
    thread_var3 = ((fft_size*thread_var2/(thread_tot*thread_var4)));
    while(thread_var3==0)
    {
      thread_var4>>=1;
      thread_var3 = ((fft_size*thread_var2/(thread_tot*thread_var4)));
    }
    
    NN=fft_size;
    thread_var1=logN;
    while(NN>thread_var2)
    {
      NN>>=1;
      thread_var1--;
    }
    pins.sewStraight(&FourierTron::CT_FFT_tiptoe, thread_tot);
    thread_var2 = thread_var4;
  }
  
  if(step_depth==0)
  {
    NN=fft_size;
    thread_var1=logN;
    while(NN>thread_tot)
    {
      NN>>=1;
      thread_var1--;
    }
    pins.sewStraight(&FourierTron::CT_FFT_march, thread_tot);
  }
  else
  {
    NN = fft_size>>1;
    thread_var3 = 0;
    while(NN != 0){NN>>=1;thread_var3++;}
    NN = thread_tot>>1;
    thread_var4 = 0;
    while(NN != 0){NN>>=1;thread_var4++;}
//     pins.sewStraight(&FourierTron::CT_FFT_march_cache, thread_tot);
    
    unsigned long int local_size = fft_size/thread_tot;
    for(unsigned long int i=0;i<thread_tot;i++)
    {
      localvec[i].resize(local_size, tempval1);
    }
    pins.sewStraight(&FourierTron::CT_FFT_march_cache_local, thread_tot);
  }
}


void FourierTron::CT_FFT_swap(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  unsigned long int max = reversal[thread_var1].size();
  for(unsigned long int i=w;i<max;i+=thread_tot)
  {
    threadvec1[w].set((*tempvecpointer1)[reversal[thread_var1][i][1]]);
    (*tempvecpointer1)[reversal[thread_var1][i][1]].set((*tempvecpointer1)[reversal[thread_var1][i][0]]);
    (*tempvecpointer1)[reversal[thread_var1][i][0]].set(threadvec1[w]);
  }
}


void FourierTron::CT_FFT_step(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int chunk=1;
  unsigned long int chunk2=2;
  unsigned long int tier=0;
  unsigned long int parts = ((fft_size>>1)>>thread_var1);
  unsigned long int basestart = w*(fft_size>>thread_var1);
  unsigned long int start=basestart;
  
  unsigned long int k1=0;
  unsigned long int k2=0;
  
  while(chunk<thread_var2)
  {
    for(unsigned long int i=0;i<chunk;i++)
    {
      twiddleVal(tier, i, threadvec1[w]);
      start=basestart;
      k1 = start + i;
      k2 = k1 + chunk;
      for(unsigned long int j=0;j<parts;j++)
      {
        mul(threadvec1[w], (*tempvecpointer1)[k2], threadvec2[w]);
        sub((*tempvecpointer1)[k1], threadvec2[w], (*tempvecpointer1)[k2]);
        (*tempvecpointer1)[k1].add(threadvec2[w]);
        k1 += chunk2;
        k2 += chunk2;
      }
    }
    tier++;
    chunk = chunk2;
    chunk2<<=1;
    parts>>=1;
  }
}


//performs a fft depth-first, until the required size of the sub-fft exceeds the cache size set by the user, at which point it proceeds breadth-first
void FourierTron::CT_FFT_step_cache(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int base_size = fft_size/thread_tot;
  unsigned long int basestart = w*base_size;
  unsigned long int breadth_steps = base_size>>step_depth;
  unsigned long int iter_skip = 1<<step_depth;
  
  //thread_var3 is the number of tiers to climb, or the log_2 of thread_var2
  unsigned long int tiers_remaining = thread_var3;
  unsigned long int steps_this_iter = step_depth;
  unsigned long int base_tier = 0;
  
  //tiers represents the level of depth in the fft
  while(tiers_remaining != 0)
  {
    if(tiers_remaining < step_depth)
    {
      steps_this_iter = tiers_remaining;
      breadth_steps = base_size>>steps_this_iter;
      iter_skip = 1<<steps_this_iter;
    }
    unsigned long int start = basestart;
    unsigned long int chunk = 1<<base_tier;
    for(unsigned long int b=0;b<breadth_steps;b++)
    {
      unsigned long int tier = base_tier;
      unsigned long int chunk2 = chunk<<1;
      unsigned long int top_parts = 1;
      unsigned long int parts = iter_skip>>1;
      //steps_this_iter is the depth traversed before going breadth-wise
      for(unsigned long int t=0;t<steps_this_iter;t++)
      {
        unsigned long int k1=0;
        unsigned long int k2=0;
        //the next two loops go breadth-wise.  The number of distinct elements swapped are meant to fit into cache
        for(unsigned long int i=0;i<top_parts;i++)
        {
          k1 = start + i*(1<<base_tier);
          k2 = k1 + chunk;
          twiddleVal(tier, k1%chunk, threadvec1[w]);
          for(unsigned long int j=0;j<parts;j++)
          {
            //finally we perform the butterfly
            mul(threadvec1[w], (*tempvecpointer1)[k2], threadvec2[w]);
            sub((*tempvecpointer1)[k1], threadvec2[w], (*tempvecpointer1)[k2]);
            (*tempvecpointer1)[k1].add(threadvec2[w]);
            k1 += chunk2;
            k2 += chunk2;
          }
        }
        tier++;
        top_parts<<=1;
        chunk = chunk2;
        chunk2<<=1;
        parts>>=1;
      }
      chunk = 1<<base_tier;
      //the logic below ensures that we run over all pairs exactly once per tier
      if(base_tier==0){start += iter_skip;}
      else if(((start+1)%chunk)==0){start += (1 + (iter_skip-1)*chunk);}
      else{start += 1;}
    }
    tiers_remaining -= steps_this_iter;
    base_tier += steps_this_iter;
  }
}


void FourierTron::CT_FFT_march_cache(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int base_size = fft_size/thread_tot;
  unsigned long int basestart = w*base_size;
  unsigned long int breadth_steps = base_size>>step_depth;
  unsigned long int iter_skip = 1<<step_depth;
  
  //thread_var3 is the number of tiers to climb
  //thread_var4 is the starting tier, or log_2 of thread_tot
  unsigned long int base_tier = thread_var4;
  unsigned long int tiers_remaining = thread_var3 - base_tier;
  unsigned long int steps_this_iter = step_depth;
  
  //tiers represents the level of depth in the fft
  while(tiers_remaining != 0)
  {
    if(tiers_remaining < step_depth)
    {
      steps_this_iter = tiers_remaining;
      breadth_steps = base_size>>steps_this_iter;
      iter_skip = 1<<steps_this_iter;
    }
    unsigned long int start = basestart;
    unsigned long int chunk = 1<<base_tier;
    for(unsigned long int b=0;b<breadth_steps;b++)
    {
      unsigned long int tier = base_tier;
      unsigned long int chunk2 = chunk<<1;
      unsigned long int top_parts = 1;
      unsigned long int parts = iter_skip>>1;
      //steps_this_iter is the depth traversed before going breadth-wise
      for(unsigned long int t=0;t<steps_this_iter;t++)
      {
        unsigned long int k1=0;
        unsigned long int k2=0;
        //the next two loops go breadth-wise.  The number of distinct elements swapped are meant to fit into cache
        for(unsigned long int i=0;i<top_parts;i++)
        {
          k1 = w + (start - basestart + i*(1<<(base_tier-thread_var4)))*thread_tot;
          k2 = k1 + chunk;
          twiddleVal(tier, k1%chunk, threadvec1[w]);
          for(unsigned long int j=0;j<parts;j++)
          {
            //finally we perform the butterfly
            mul(threadvec1[w], (*tempvecpointer1)[k2], threadvec2[w]);
            sub((*tempvecpointer1)[k1], threadvec2[w], (*tempvecpointer1)[k2]);
            (*tempvecpointer1)[k1].add(threadvec2[w]);
            k1 += chunk2;
            k2 += chunk2;
          }
        }
        tier++;
        top_parts<<=1;
        chunk = chunk2;
        chunk2<<=1;
        parts>>=1;
      }
      chunk = 1<<base_tier;
      //the logic below ensures that we run over all pairs exactly once per tier
      if(((start+1)%(chunk>>thread_var4))==0){start += (1 + (iter_skip-1)*(chunk>>thread_var4));}
      else{start += 1;}
    }
    tiers_remaining -= steps_this_iter;
    base_tier += steps_this_iter;
  }
}


void FourierTron::CT_FFT_march_cache_local(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  //initialize localvec
  unsigned long int offset = 0;
  for(unsigned long int i=0;i<localvec[w].size();i++)
  {
    localvec[w][i].set((*tempvecpointer1)[w + offset]);
    offset += thread_tot;
  }
  unsigned long int base_size = fft_size/thread_tot;
  unsigned long int basestart = w*base_size;
  unsigned long int breadth_steps = base_size>>step_depth;
  unsigned long int iter_skip = 1<<step_depth;
  
  //thread_var3 is the number of tiers to climb
  //thread_var4 is the starting tier, or log_2 of thread_tot
  unsigned long int base_tier = thread_var4;
  unsigned long int tiers_remaining = thread_var3 - base_tier;
  unsigned long int steps_this_iter = step_depth;
  
  //tiers represents the level of depth in the fft
  while(tiers_remaining != 0)
  {
    if(tiers_remaining < step_depth)
    {
      steps_this_iter = tiers_remaining;
      breadth_steps = base_size>>steps_this_iter;
      iter_skip = 1<<steps_this_iter;
    }
    unsigned long int start = basestart;
    unsigned long int chunk = (1<<base_tier)/thread_tot;
    for(unsigned long int b=0;b<breadth_steps;b++)
    {
      unsigned long int tier = base_tier;
      unsigned long int chunk2 = chunk<<1;
      unsigned long int top_parts = 1;
      unsigned long int parts = iter_skip>>1;
      //steps_this_iter is the depth traversed before going breadth-wise
      for(unsigned long int t=0;t<steps_this_iter;t++)
      {
        unsigned long int k1=0;
        unsigned long int k2=0;
        //the next two loops go breadth-wise.  The number of distinct elements swapped are meant to fit into cache
        for(unsigned long int i=0;i<top_parts;i++)
        {
          //w + n*thread_tot --> n
          //chunk --> chunk/thread_tot
//           k1 = w + (start - basestart + i*(1<<(base_tier-thread_var4)))*thread_tot;
          k1 = (start - basestart + i*(1<<(base_tier-thread_var4)));
          k2 = k1 + chunk;
          twiddleVal(tier, (k1*thread_tot + w)%(chunk*thread_tot), threadvec1[w]);
          for(unsigned long int j=0;j<parts;j++)
          {
            //finally we perform the butterfly
            mul(threadvec1[w], localvec[w][k2], threadvec2[w]);
            sub(localvec[w][k1], threadvec2[w], localvec[w][k2]);
            localvec[w][k1].add(threadvec2[w]);
            k1 += chunk2;
            k2 += chunk2;
          }
        }
        tier++;
        top_parts<<=1;
        chunk = chunk2;
        chunk2<<=1;
        parts>>=1;
      }
      chunk = (1<<base_tier);
      //the logic below ensures that we run over all pairs exactly once per tier
      if(((start+1)%(chunk>>thread_var4))==0){start += (1 + (iter_skip-1)*(chunk>>thread_var4));}
      else{start += 1;}
      chunk/=thread_tot;
    }
    tiers_remaining -= steps_this_iter;
    base_tier += steps_this_iter;
  }
  offset=0;
  for(unsigned long int i=0;i<localvec[w].size();i++)
  {
    (*tempvecpointer1)[w + offset].set(localvec[w][i]);
    offset += thread_tot;
  }
}


void FourierTron::CT_FFT_tiptoe(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  unsigned long int w1 = w/thread_var2;
  unsigned long int w2 = w%thread_var2;
  unsigned long int W = w1*thread_var4 + w2;
  
  unsigned long int jump=thread_var2;
  unsigned long int tier=thread_var1;
  unsigned long int i2=1;
  unsigned long int i1 = (thread_var4/thread_var2)>>1;
  unsigned long int skip = thread_var2<<1;
  
  unsigned long int k1=0;
  unsigned long int k2=0;
  unsigned long int k1_j,k1_i;
  
  unsigned long int niter = thread_var3;
  unsigned long int nskip = fft_size/niter;
  
  while(jump<thread_var4)
  {
    k1_j = W;
    for(unsigned long int j=0;j<i2;j++)
    {
      twiddleVal(tier, (W + j*thread_var2)%(jump<<1), threadvec1[w]);
      k1_i = 0;
      for(unsigned long int i=0;i<i1;i++)
      {
        k1 = k1_j + k1_i;
        for(unsigned long int k=0;k<niter;k++)
        {
          k2 = k1 + jump;
          
          mul(threadvec1[w], (*tempvecpointer1)[k2], threadvec2[w]);
          sub((*tempvecpointer1)[k1], threadvec2[w], (*tempvecpointer1)[k2]);
          (*tempvecpointer1)[k1].add(threadvec2[w]);
          
          k1 += nskip;
        }
        k1_i += skip;
      }
      k1_j += thread_var2;
    }
    jump<<=1;
    skip<<=1;
    tier++;
    i2<<=1;
    i1>>=1;
  }
}


void FourierTron::CT_FFT_march(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  unsigned long int ntemp = thread_tot;
  unsigned long int ntemp2 = (ntemp<<1);
  unsigned long int part = ((fft_size>>(thread_var1))>>1);
  unsigned long int start=0;
  unsigned long int tier = thread_var1;
  
  unsigned long int k1=0;
  unsigned long int k2=0;
  
  while(part!=0)
  {
    for(unsigned long int i=w;i<ntemp;i+=thread_tot)
    {
      twiddleVal(tier, i, threadvec1[w]);
      start=0;
      k1 = start + i;
      k2 = k1 + ntemp;
      for(unsigned long int j=0;j<part;j++)
      {
        mul(threadvec1[w], (*tempvecpointer1)[k2], threadvec2[w]);
        sub((*tempvecpointer1)[k1], threadvec2[w], (*tempvecpointer1)[k2]);
        (*tempvecpointer1)[k1].add(threadvec2[w]);
        
        k1 += ntemp2;
        k2 += ntemp2;
      }
    }
    tier++;
    ntemp=ntemp2;
    ntemp2<<=1;
    part>>=1;
  }
}


void FourierTron::CT_FFT(const vector<cxfloat> &input, vector<cxfloat> &output, unsigned long int N)
{
  //check that input.size() >= N
  try{if(input.size() < N){throw (char*)("input vector size smaller than third argument");}}
  catch(char* str){cout<<"Exception from FourierTron::CT_FFT : "<<str<<endl;throw;return;}
  
  if(output.size() < N){output.resize(N, tempval1);}
  for(unsigned long int i=0;i<N;i++){output[i].set(input[i]);}
  CT_FFT(output, N);
}














